##############################################################
#Figure 1: Spatial Distribution of Hurricane Paths 2003-2014
##############################################################


## Figure 1 (a) 

fips <- read.csv("https://www2.census.gov/geo/docs/reference/codes/files/national_county.txt", header=F, sep=",", 
                 colClasses = c(rep("character",5)))

fips$FIPS_initial<- with(fips, paste(V2, V3))
fips$FIPS <- gsub(" ", "", fips$FIPS, fixed = TRUE)

fips$V1[fips$V1=="AL"]<- "Alabama"	            
fips$V1[fips$V1=="AK"]<- "Alaska"	            
fips$V1[fips$V1=="AZ"]<- "Arizona"	            
fips$V1[fips$V1=="AR"]<- "Arkansas"	        
fips$V1[fips$V1=="CA"]<- "California"	        
fips$V1[fips$V1=="CO"]<- "Colorado"	        
fips$V1[fips$V1=="CT"]<- "Connecticut"	        
fips$V1[fips$V1=="DE"]<- "Delaware"	        
fips$V1[fips$V1=="FL"]<- "Florida"	            
fips$V1[fips$V1=="GA"]<- "Georgia"	            
fips$V1[fips$V1=="HI"]<- "Hawaii"	            
fips$V1[fips$V1=="ID"]<- "Idaho"	            
fips$V1[fips$V1=="IL"]<- "Illinois"	        
fips$V1[fips$V1=="IN"]<- "Indiana"	            
fips$V1[fips$V1=="IA"]<- "Iowa"	            
fips$V1[fips$V1=="KS"]<- "Kansas"	            
fips$V1[fips$V1=="KY"]<- "Kentucky"	        
fips$V1[fips$V1=="LA"]<- "Louisiana"	        
fips$V1[fips$V1=="ME"]<- "Maine"	            
fips$V1[fips$V1=="MD"]<- "Maryland"	        
fips$V1[fips$V1=="MA"]<- "Massachusetts"	    
fips$V1[fips$V1=="MI"]<- "Michigan"	        
fips$V1[fips$V1=="MN"]<- "Minnesota"	        
fips$V1[fips$V1=="MS"]<- "Mississippi"	        
fips$V1[fips$V1=="MO"]<- "Missouri"	        
fips$V1[fips$V1=="MT"]<- "Montana"	            
fips$V1[fips$V1=="NE"]<- "Nebraska"	   		
fips$V1[fips$V1=="NV"]<- "Nevada"				
fips$V1[fips$V1=="NH"]<- "New Hampshire"		
fips$V1[fips$V1=="NJ"]<- "New Jersey"			
fips$V1[fips$V1=="NM"]<- "New Mexico"			
fips$V1[fips$V1=="NY"]<- "New York"			
fips$V1[fips$V1=="NC"]<- "North Carolina"		
fips$V1[fips$V1=="ND"]<- "North Dakota"		
fips$V1[fips$V1=="OH"]<- "Ohio"				
fips$V1[fips$V1=="OK"]<- "Oklahoma"			
fips$V1[fips$V1=="OR"]<- "Oregon"				
fips$V1[fips$V1=="PA"]<- "Pennsylvania"		
fips$V1[fips$V1=="RI"]<- "Rhode Island"		
fips$V1[fips$V1=="SC"]<- "South Carolina"		
fips$V1[fips$V1=="SD"]<- "South Dakota"		
fips$V1[fips$V1=="TN"]<- "Tennessee"			
fips$V1[fips$V1=="TX"]<- "Texas"				
fips$V1[fips$V1=="UT"]<- "Utah"				
fips$V1[fips$V1=="VT"]<- "Vermont"				
fips$V1[fips$V1=="VA"]<- "Virginia"			
fips$V1[fips$V1=="WA"]<- "Washington"			
fips$V1[fips$V1=="WV"]<- "West Virginia"		
fips$V1[fips$V1=="WI"]<- "Wisconsin"			
fips$V1[fips$V1=="WY"]<- "Wyoming"				

fips$state <- sapply(fips$V1, tolower)

fips$county1 <- fips$V4

fips$county1 = substr(fips$county1,1,nchar(fips$county1)-6)

fips$county2 <- sapply(fips$county1, tolower)
fips$county <- trimws(fips$county2, "r")
head(fips)
new_fips <- fips[c(-1,-2, -3, -4, -5, -6, -9, -10)]

new_fips$countystate<- with(new_fips, paste(state, ",", county))
new_fips$countystate <- gsub(" ", "", new_fips$countystate, fixed = TRUE)
new_fips$county <- new_fips$countystate
head(new_fips)
fips_data <- new_fips[c(-2,-3)]
head(fips_data)



county <- read.dta13("./input/county_final.dta")
head(county)


fips_final = merge(fips_data, county, by.x="countystate", by.y="countystate", all =TRUE)

hist(fips_final$no_hur_exposed)

fips_final$break2 <- cut(fips_final$no_hur_exposed, 
                         breaks=c(0, 2, 4, 6, 8, 11), 
                         include.lowest=TRUE)



world <- ne_countries(scale='medium',returnclass = 'sf')
us <- ggcounty.us() 



ggplot(data = world) +
  geom_map(data=fips_final[fips_final$distance<2500,], map=us$map, aes(map_id=FIPS, fill=distance),  size=0.125) + 
  geom_polygon(data = map_data("state"), aes(x = long, y = lat, group = group), 
               fill = "antiquewhite1", colour = "white", size = 0.3, alpha=.5) + 
  annotate(geom = "text", x = -90.5, y = 27, label = "Gulf of Mexico", 
           color = "grey22", size = 4.5) +
  coord_sf(xlim = c(-105, -74), ylim = c(25, 40)) +
  xlab("Longitude")+ ylab("Latitude")+
  theme(panel.grid.major = element_line(colour = gray(0.5), linetype = "dashed", size = 0.5),
        panel.background = element_rect(fill = "aliceblue"),
        panel.border = element_rect(fill = "NA"))  + 
  scale_fill_gradientn(colours = c("grey80","red", "green", "cyan", "antiquewhite1"),
                       values = scales::rescale(c(0, 30, 200, 500, 1000, 2500)), name="Distance")  

ggsave(file="./figures/fig1.pdf")



### Figure 1 (b)


myhurricane <- read.csv("~/Dropbox/Replication/JoP/input/myhurricanes.csv")

myhurricane$ID = as.factor(paste(myhurricane$Name, myhurricane$Season, 
                                 sep = "."))
myhurricane$Name = as.factor(myhurricane$Name)


#Tropical Depression

myhurricane$Category[myhurricane$Wind.WMO.<=33] <- -2 

#Tropical Storm

myhurricane$Category[myhurricane$Wind.WMO.>=34 & myhurricane$Wind.WMO.<=63] <- -1 

myhurricane$Category[myhurricane$Wind.WMO.>=64 & myhurricane$Wind.WMO.<=82] <- 1 
myhurricane$Category[myhurricane$Wind.WMO.>=83 & myhurricane$Wind.WMO.<=95] <- 2
myhurricane$Category[myhurricane$Wind.WMO.>=96 & myhurricane$Wind.WMO.<=112] <- 3
myhurricane$Category[myhurricane$Wind.WMO.>=113 & myhurricane$Wind.WMO.<=136] <- 4
myhurricane$Category[myhurricane$Wind.WMO.>=137] <- 5
#windwmo measured in kt 


ggplot(myhurricane, aes(x = Longitude, y = Latitude, group = ID)) +
  geom_polygon(data = map_data("state"), aes(x = long, y = lat, group = group), 
               fill = "antiquewhite1", colour = "grey", size = 0.2, alpha=.3) + 
  geom_path(data = myhurricane, aes(group = ID, colour = as.factor(Category)), size=.8) +
  xlim(-125,-67) + ylim(21, 50) + 
  labs(x = "", y = "", colour = "Category") + 
  theme(panel.background = element_blank(),
        axis.line = element_blank(),
        axis.title = element_blank(),
        axis.text=element_blank(),
        axis.ticks = element_blank(),
        legend.key = element_rect(fill = "transparent", colour = "transparent"),
        legend.position = "right") + 
  scale_color_manual(values = c("skyblue", "deepskyblue", "lightgoldenrod1", "goldenrod1", "indianred2", "brown1", "brown3"), labels=c("TD", "TS", "1", "2", "3", "4", "5")) 

ggsave(file="./figures/fig1b.pdf")


